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ABSTRACT 

Loss of equilibrium of magnetic flux ropes is a leading candidate for the origin of solar coronal mass 
ejections (CMEs). The aim of this paper is to explore to what extent this mechanism can account 
for the initiation of CMEs in the global context. A simplified MHD model for the global coronal 
magnetic field evolution in response to flux emergence and shearing by large-scale surface motions 
is described and motivated. Using automated algorithms for detecting flux ropes and ejections in 
the global magnetic model, the effects of key simulation parameters on the formation of flux ropes 
and the number of ejections are considered, over a 177-day period in 1999. These key parameters 
include the magnitude and sign of magnetic helicity emerging in active regions, and coronal diffusion. 
The number of flux ropes found in the simulation at any one time fluctuates between about 28 and 
48, sustained by the emergence of new bipolar regions, but with no systematic dependence on the 
helicity of these regions. However, the emerging helicity does affect the rate of flux rope ejections, 
which doubles from 0.67 per day if the bipoles emerge untwisted to 1.28 per day in the run with 
greatest emerging twist. The number of ejections in the simulation is also increased by 20% - 30% by 
choosing the majority sign of emerging bipole helicity in each hemisphere, or by halving the turbulent 
diffusivity in the corona. For reasonable parameter choices, the model produces approximately 50% 
of the observed CME rate. This indicates that the formation and loss of equilibrium of flux ropes may 
be a key clement in explaining a significant fraction of observed CMEs. 

Subject headings: Sun: corona — Sun: coronal mass ejections — Sun: evolution — Sun: magnetic 
fields 
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1. INTRODUCTION 

Coronal mass ejections (CMEs) from the Sun are im- 
portant drivers of space weather ( |Schwenn 2006 1 . They 
are believed fundamental to the long-term evolution of 
the global solar magne tic field, both by shedding mag 



netic flux and helicity (Low 19961 Bieber & Rust 



Lynch et al.|2005| , and by their role in the corona 



1995 



2001 Owens et al. 



2007). 



mag- 



netic field reversal every 11 years (e.g., Zhang & Low 



Avariety of models have been pro posed for the CME 
initiation process (see the reviews by| Forbes|2000| |Klim-| 
|chuk| [2001 Low 200T|. Although models differ in their 
basic magnetic field configuration and in how the erup- 
tion is initiated, there is agreement that t he driving en - 
ergy must originate in the magnetic field ( |Forbes|[ 2000 ) . 
The majority of recent work o n CME initiation favours 
"storage-and-release" models ( |Klimchuk| |2001[ ), where 
free magnetic energy is built up gradually in a quasi- 
static evolution, before sudden release in a highly dy- 
namic eruptive phase. There are essentially two pre- 
emption configuratio ns invoked to store the free energy: 
sheared arcades (e.g., |Antiochos et al.|1999 ) or magnetic 
flux ropes. The latter are twisted structures of helical 
magnetic fields, anchored at both ends in the photo- 
sphere. Note that the degree of twist may be low, and 
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flux ropes with less than one turn are essentially the same 
as sheared arcades. Indeed, many C ME models wor k 
with either a flux rope or an arcade (Klimchuk 20011); 



what matters is that sufficient magnetic pressure builds 
in the structure to overcome the stabilizing magnetic ten- 
sion of the overlying field. There is, however, a consider- 
able body of evidence that at least a third of CMEs have 
a flux rope structure, ranging fro m white-light corona- 



graph and EUV observations (e.g., Illing fc Hundhausen 



1986 Dere et al. 1999 Gibson et al. 2006) to in site obser- 
vations of rotating ma gnetic field patterns in i ntcrplanc- 
tary magnetic clouds (Klein & Burlaga 19821. Clearly 



twisted structures are observed in i mages of erupting 
prominences at various wavelengths (Tandbcrg-Hanssen 



1974| |Plunkett et al-pOOO] |Green et al ||2007p . It is pos 
sible that flux r opes observed in CM Es may be formed 



during eruption (Gosling et al. 19951 rather than existing 
previously in a quiescent state. However, the existence 
of coronal cavities before lift-off is s trong evidence for a 
pre-existing flux rope in many cases ( Gibson et al.|2006[ ), 
because the flux rope can explain both the lower density 
of the cavity (due to strong magnetic pres sure) , and the 



suppo rt of filament mass inside the cavity ( Gibson & Fan 
[2006| . 

Previous models for CME initiation have considered 
only a single event in a simplified magnetic configura- 
tion, in order to study the basic physical process. In 
this paper, however, we begin to address the fundamen- 
tal question of where and when CMEs occur in the global 
context. A quasi-static model is used to study the build 



2 



Yeates et al. 



up of axial flux in magnetic flux ropes, using observed 
photospheric flux distributions as input to drive simu- 
lations of the global coronal magnetic field. In essence, 
this model is an extension of the two-dimensional (2D) 



catast rophe models considered by|van Tend fc Kuperus 



([1978 1 



et al 



and many authors since (see the review by Lin 
2003). These 2D catastrophe models describe the 



equilibrium curve of a flux rope as some control param- 
eter, typically describing the background magnetic field, 
is varied. Eventually, a nose point in the equilibrium 
curve may be reached beyond which no equilibrium is 
possible. In the present 3D global model, many flux 
ropes form at different locations on the Sun, and may 
or may not lose equilibrium as they evolve. The location 
of each flux rope and evolution of its background field are 
determined self-consistently based on the observed pho- 
tospheric magnetic field, rather than being prescribed as 
in the 2D models. Therefore, in this paper, the evolution 
of the coronal field leading to the formation and ejection 
of flux ropes is constrained by observations. 

The global model described in this paper differs from 
other existing models of the global magnetic field in 
the solar corona. The majority of such models have 
used potential-field extrapolations with the radial mag- 
netic field in the solar photosphere as a lower boundary 
condition, and an upper source surface at r = 2 .5-R^ 



study the time evolution of the corona with potential- 
field models, a sequ ence of independent extrapolations 



has been used (e.g., [Wa ng et al. 2002] |Mackay fc Lock 
wood 20021 |Schrijver & Derosa 2003 ) . Unfortunately, the 



potential field suffers from the fundamental limitation 
of having the minimum magnetic energy for the given 
boundary conditions, thus containing no free magnetic 
energy to drive CMEs. More general extrapolation meth- 
ods are under development for the global corona, includ- 
ing nonlinear f orce-free fields base d on full-disk vector 
magnetogr ams ( Wi egclmann 20 071) and f ull-MHD solu- 



tions (e.g., |Mikic et al.||1999| |Riley et al ||2006| ). How- 
ever, such models are too computationally intensive to 
simulate the long-term evolution (over periods of months 
or longer), so we have developed a new technique using 
a coup led surface flux transport and magneto-frictional 
model (lYeates et al. 120071 |2008al hereafter "Paper I" and 



"Paper II"). This allows us to insert active regions with 
non-zero magnetic helicity into the global corona, and to 
follow the evolution of twisted and sheared non-potential 
magnetic fields through a sequence of equilibria as they 
are driven by photospheric motions. Like the 2D catas- 
trophe models, a quasi-static approximation of the evo- 
lution is used, so that the model can follow the build up 
of flux ropes and their loss of equilibrium, but not the 
subsequent dynamics of the eruption, which would re- 
quire full time-dependent MHD simulations. Note that 
coronal flux ropes form not only within active regions, 
but also between regio ns, consistent with the observa - 



tions of solar filaments (Tang|1987 Mackay et al 



2008) 



orma- 



The research presented here, m considering the 
tion and ejection of flux ropes, has been motivated by 
the success of the model in reproducing the hemispheric 
chirality pattern of sheared magnetic fields in filaments 
(Paper II) as non-potential fields of solar filaments are 
known to be related to CMEs. 
The aim of this paper is to describe the model and its 



( |Altschuler fc Newkirk|[l969| |Schatten et al ||1969[ ). To h en e t al 



physical motivation, and to demonstrate how the forma- 
tion and evolution of magnetic flux ropes in the model 
depends on certain parameters, in particular the mag- 
netic helicity of emerging regions. We physically test 
whether flux cancellation, as simulated in this model, 
can produce enough flux rope ejections to match the ob- 
served CME rate, or whether other physical mechanisms 
are required. A detailed comparison of our simulation 
with the observed locations of CMEs will be the subject 
of a future paper. The organization of this paper is as 
follows. The model and its physical assumptions are de- 
scribed in Section [2j In Section [3] we describe the basic 
process of flux rope formation and ejection that occurs 
many times in the global model. Automated techniques 
for detecting flux ropes and their ejection in the global 
simulations are described in Section [4] then applied in 
Section [5] to consider the effect of simulation parameters 
on the flux ropes formed and on the rate of ejections. 
In Section [6] we summarize our conclusions, and suggest 
future extensions of the work. 

2. CORONAL EVOLUTION MODEL 

A full treatment of the dynamical evolution of the 
global solar corona over weeks to months is computa- 
tionally prohibitive. Instead, we use the coupled flux 
transport a nd ma gneto-frictional model of van Ballcgooi- 



(2000), which describes the evolution of the 
large-scale mean magnetic field under certain assump- 
tions. The basic physics of the model has been studied in 
a series of papers with the aim of understanding the chi- 
rality of observed solar filaments. These papers include 
detailed parameter st udies of a pair of interacting bipo 



lar magnetic r egions (Mackay & van Ballegooijen 2001 
in addition to simulations drive n by ob 



|2005| 2006a|b[ 

served pho tospheric magnetic fields, both l ocal ([Mackay 
let al.|[2000[ ), and global (Paper I, Paper II, |Yeates et al. 
2008b), as in this paper. In this section we summarize 
the basic assumptions of this model, before describing 
our particular numerical implementation. We describe 
the observational input required for the simulations to 
accurately reproduce the photospheric and coronal mag- 
netic fields. 

2.1. Basic Physical Assumptions 

To model the evolution of the large-scale magnetic field 
in the solar corona, we make the following simplifying 
assumptions: 

1. The magnetic field in the solar corona may be 
separated into a large-scale, mean, component on 
lengthscales ~ 50 Mm, and a fluctuating compo- 
nent on lengthscales ~ 1 Mm. 

2. On the timescale of large-scale photospheric mo- 
tions, the large-scale coronal magnetic field is in 
equilibrium, except perhaps during flares or CME 
eruptions. 

3. Above the photosphere, and in the low corona up 
to 2.5i?0, magnetic forces dominate over plasma 
forces (i.e., it is a low- f3 plasma). In the photo- 
sphere magnetic flux tubes are passively advected. 

4. The large-scale magnetic field in the photosphere 
is predominantly radial, and its horizontal length- 
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scale is larger than that of the supergranular con- 
vection. 

5. For purposes of the large-scale magnetic field evo- 
lution, emerging active regions may be represented 
by twisted magnetic bipoles, with an idealized 
form. 

Assumption 1 allow s the application of mean field the- 
ory to t he corona (Seehafer 1994 |van Ballegooijen et 



al.||2000 [). In that case, only the mean magnetic field 
is solved for, and the effect of significant small-scale 
structure — such as braiding and current sheets produced 
by interaction with convective flows in the photosphere 
is included through a turbulent electro- 



( Parker 


19721 


motive 


orce ( 



orce (e.m.f.). The mean-field induction equation 
takes the form 

dA 



dt 



v x B + £ , 



(1) 



where Ao(r, t) is the vector potential for the mean mag- 
netic field Bo = V x Ao, vo(r,t) is the mean plasma 
velocity, and £o(r,t) is the turbulent e.m.f.. 

Assumption 2 is valid because magnetic disturbances 
from the photosphere propagate up into the corona 
at Alfven speeds o f the order 1000 km s -1 or greater 



(Regnier et al. 2008| , well above the maximum flow speed 
of surface rotation (of the order 0.2kms _1 ). The evolu- 
tion of the coronal magnetic field may then be under- 
stood as a continual distortion of existing equilibrium by 
footpoint motions, with s ubsequent fast r elaxation to a 
neighbouring equilibrium ( |Seehafer]|1994[ ) . By assump- 
tion 3, the form of this equilibrium is well approximated 
by a force-free field j x B = 0, where jo = V x B /Vo 
is the mean current density. 

In contrast to the corona, the photospheric field is as- 
sumed to be passively advected by plasma flows. As- 
sumption 4 then justifies use of the standard surface flux 
transport model for the evoluti on of the radia l compo- 
nent Bo r on the solar surface (Sheeley 2005 and ref- 
erences therein). That the large-scale photospheric field 
should be approximately rad ial is supported both by vec - 
tor magnetic measurements ( Martinez Pil let et al.|1997 1 , 
and by theoretical considerations, due to the concentra- 
tion of convection zone magnetic flux into iso lated kilo- 
gauss flux tubes (van Ballegooijen & Mackay 2007). In 
the surface flux transport model, B 0r is advected by the 
large-scale motions of differential rotation and meridional 
circulation, and the influence of supergranular convection 
on thi s large-scale fie ld is modelled by a surface diffusion 
term ( Leighton|1964 1. In the model described in this pa- 
per, surface fiux transport, coupled with the emergence 
of new active regions, acts as the lower boundary condi- 
tion to drive the coronal field evolution. 

Assumption 5 is the most uncertain, due to the ide- 
alized nature of the magnetic bipoles. Previous studies 
have shown that such idealized bipoles are adequate for 
reproducing the distribution of large-scale radial field Bq t . 
on the solar surface (Paper I). The main uncertainty is 
the inclusion of twist. While our implementation in this 
paper is simplistic, it is better than not including any 
twist. Such twist has been shown to be required both 
by vector magnetic field measurements ( |Demoulin|2007| ) 
and by our previous simulations (Paper 11). Part of the 
aim of the present work is to determine whether such a 



simplified model can yield useful information about the 
global magnetic field, or whether more detailed modelling 
of individual regions is necessary to capture even the sim- 
plest large-scale phenomena. 

2.2. Numerical Implementation 

As in Paper II, we solve Equation Q in a domain 
extending from 0° to 360° in longitude, -80° to 80° in 
latitude, and Rq to 2.5R Q in radius. We assume the 
turbulent e.m.f. to take the form 



£o = -r)ja, 
with the turbulent diffusivity 

V = Vo(l + 0.2 



I Jo 



|Bo| 



(2) 



(3) 



as in Mackay & van Ballegooijen (2006a I. The first 
term is a uniform background value and the second term 
is an enhancement in regions of strong current density. 
This e nhancement was introduced becau se earlier simula- 
tions (Mackay & van Ballcgooije n|2001 1 produced highly 
twisted^ structures which are t ypically not observed (Bo- 
bra et al. |[2008] |Su et al j]2009| ). 

To follow the evolution through a sequence of nonlinear 
force-free states i n res ponse to flux emergence and surface 



shearing (Section 2.1 1, the MHD momentum equat ion is 
approximat ed by the magneto-frictional method (Yang 
et al.||198"6| , setting 



1 Jo x B , 
v " = - B 2 + v out(r)r- 



(4) 



The second term in Equation Q is a radial outflow im- 
posed near the upper boundary, where it simulates the 
effect of the so lar wind in opening up field lines in the ra- 
dial direction ( Mackay fc van Ballegooijen|2006a ) . It has 
a peak value vq = 100kms _i and falls off exponentially 
away from the upper boundary, so that only structures 
near this boundary are affected. 

We assume closed boundaries in the latitudinal direc- 
tion and an open boundary at the upper source surface 
(r = 2.5Rq), where the radial outflow ensures that the 
magnetic field is radial. At the lower boundary (r = Rq) 
the radial magnetic field Bq t is evolved according to the 
surface flux transport model with flux emergence, as de- 
scribed in Paper I. 

In this paper, we simulate the solar corona over 177 
days in 1999, the same period as Paper II. The initial con- 
dition is a potential field extrapolation for 1999 April 16 
(day of year 106) taken from a Kitt Peak synoptic magne- 
togram corrected for differential rotation (Paper I). The 
coronal magnetic field is then evolved continuously using 
Equations (HI and Q , with surface flux transport on the 
lower boundary ana flux emergence, until 1999 October 
10 (day of year 283). During this period, 119 new active 
regions are inserted in the form of 3D twisted magnetic 
bipoles, to be described in Section [273] The model then 
follows the quasi-static build-up of magnetic energy in 
the solar corona. 

2.3. Observational Input 

The emerging magnetic bipoles, which constitute the 
primary observational input to our simulations, are 
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Fig. 1. — Comparison between observed synoptic magnetograms (panels a, b) and simulated radial magnetic field Bq t on the photosphere 
(panels c, d), on day 147 and day 283. White shading indicates positive flux and black shading negative, with a saturation level of 
±25 G. The observed magnetograms (a, b) are corrected for differential rotation and smoothed for easier comparison. On the simulated 
magnetograms (c, d), the zero contours from the observed magnetograms (a, b) are overlayed in white. To the right of the white dashed 
line, new regions from the next rotation will appear in (c) but not in (a), because they are inserted 7 days before observation at central 
meridian. 



based on active regions observed in synoptic normal- 
component magnetograms from the National Solar Ob- 
servatory (NSO) at Kitt Peak. For each of 119 bipolcs 
emerged during this 177-day period, we recorded the 
day of central meridian passage, latitude, size, tilt an- 
gle, and magnetic flux from the Kitt Peak data (Pa- 
per I). Rather than using the actual magnetic field dis- 
tribution of each region, we insert idealized magnetic 
bipoles — of the mathematical form given in Paper II — 
into the simulation, with parameters chosen to match 
the above observed properties. The advantage of this 
approach is that the bipoles may readily be inserted in 
3D with photospheric and coronal components, and non- 
zero magnetic helicity. To model each observed active 
region more accurately would require the availability of 
vector magnetograms, and an extrapolation to model the 
unknown coronal magnetic field. Methods are under de- 
velopment to extrapolate nonlinear force-free fields from 
photospheric vector magnetograms, but there are still 



seriou s problems with this procedure (see DeRosa et al. 
20091, apart from the major computational cost for a 
large number of regions. Our idealized bipoles cannot 
reproduce the detailed structure and dynamics within 
individual active regions. However, the model has been 
shown to reproduce both the long-term global disper- 
sal of active regions across the solar surface (Paper I) 
and the statistical hemispheric pattern of filament chiral- 
ity (Paper II). The reproduction of the observed surface 
magnetic field is illustrated in Figure fl] which shows the 
observed and simulated Bq t (Rq, 8, cp)i on day 147 (after 
41 days' evolution) and on day 283 (at the end of the sim- 
ulation) . The accuracyof the simulated field may be seen 
from Figures [TJc) andfljd), where the zero contours from 
the observed magnetograms are overlayed on the simu- 



lated flux distribution. We now consider two aspects of 
the idealized bipoles which are poorly determined from 
observations. 

2.3.1. Bipole Insertion 

Because we only observe one side of the Sun, the exact 
date of emergence of most active regions is unknown, as 
are the details of the flux emergence process. For simplic- 
ity, we insert each region instantaneously, on an arbitrary 
day which we choose to be 7 days before central-meridian 
passage of the region. As detailed in Paper I, the prop- 
erties of the bipole are "evolved" back in time so as to 
reproduce the observed properties when passing central 
meridian. Rather than simply adding the bipole vector 
potential to any pre-existing magnetic field, we "sweep" 
away pre-existing field from the insertion region (Fig. 3 
of Paper II) . This creates a more realistic end-state for 
the emerged region as an independent flux system, and 
prevents the formation of disconnected flux in the corona. 
It also crudely models the distortion of pre-existing coro- 
nal field by a newly-emerging flux t ube, as found in MHD 
simulations of flux emergence (e.ff., |Yokoyama fc Shibata 
|1996| |Krall et aL] [T998 ) . Of course, the sweeping proce- 
dure has an impact on the structure of sheared magnetic 
fields around the bipole, and thus on some of the mag- 
netic flux ropes considered in this paper. However, in 
a detailed study of the origin of sheared fields at 109 



filam ent locations in the simulation (|Yeates fc Mackay 



2009 hereafter "Paper III"), the sweeping was found to 



be important in only 19 cases. 



2.3.2. 



Bipole Twist 

The 3D magnetic bipoles are given a non-zero magnetic 
helicity through a parameter [3 (defined in Equations 6 
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TABLE 1 

Summary of parameters in different simulation runs. 



Run 


P in N. hemisphere 


P in S. hemisphere 


r) (km 2 s 1 ) 


dq (kms 1 ) 


AN 


No emerging regions 


45 


100 


Am6 


0.6 


-0.6 


45 


100 


Am4 


0.4 


-0.4 


45 


100 


Am2 


0.2 


-0.2 


45 


100 


AO 








45 


100 


A2 


-0.2 


0.2 


45 


100 


Ad 


-0.4 


0.1 


45 


100 


A6 


-0.6 


0.6 


45 


100 


D4 


-0.4 


0.4 


22.5 


100 


VI 


-0.4 


0.1 


45 


50 




Fig. 2. — Effect of the twist parameter /3 on the shape of a 
3D magnetic bipole with (a) j3 = 0.2 and (b) /3 = —0.6. Grey 
shading and thin contours show the radial magnetic field Brjr on 
the photosphere (white/solid contours for positive, grey/dashed 
contours for negative) . Thick lines show selected coronal magnetic 
field lines. 



at sites of flux emergence, but also between the bipoles 
where flux was cancelling. The magnitude of (3 affected 
the formation rate of magnetic flux ropes, and hence the 
time until loss of equilibrium, but did not affect whether 
or not a flux rope formed. Ma ckay fc van Ballegooijen] 
( 2005 1 showed that the sign of (3 also affects the chiral- 



ity type (direction of shear) in this simple configuration. 
Subsequent studies of the chirality generated in global 
simulations have found that both emerging helicity and 
shearing by surface motions can be important at differ- 
ent locations on the Sun (Paper III) . For the 109 filament 
locations examined in that study, only 32% changed chi- 
rality when the sign of (3 in nearby regions was reversed. 

Ideally, the magnitude of helicity in each bipole would 
be based on vector magnetograph data for the real active 
region in 3D space, but such measurements are unavail- 
able in the corona. Although a number of authors have 
attempted to estimate helicity in active regions using vec- 
tor magnetic fie ld measurements in the photosphere (see 



Demoulin 2007), the horizontal magnetic field compo- 
nents have high uncertainties — including a 180° degree 
ambiguity ( jMetcalf et al.||2006[ ) — which are exacerbated 
by taking derivatives to compute the vertical current den- 
sity. Thus we are unable to determine observationally 
the best value of (3 to model each active region. We 
therefore assume that all bipoles in each hemisphere have 
the same value of (3, and run simulations with different 
values to determine the effect this has on the behav- 
ior of the model. Note that taking the same value of 
(3 for all regions in each hemisphere does not lead to 
uniform magnetic helicity in space as the bipoles have 
di fferent sizes an d field s trengths. This is well illustrated 
(2008bl in plots of the "current hclic- 
Bq/Bq, a measure of the local twist 
One seemingly robust 



by Yeates et al 
V 



ity* a 



B 



X JDQ 

of field lines in a force-free field 
feature of magnetic helicity in the solar corona is the 
hemispheric pattern, whereby the majority of active re- 
gions in the Northern hemisphere have negative helicity, 
and in the Southern hemisphere, positive helicity. This 



and 8 of Paper II) , which describes the initial twist of the 
field lines in the corona, and does not affect the bipole's 
radial magnetic footprint on the solar surface. Figure [2] 
shows the structure of a bipole with two different mag- 
nitudes and signs of (3, taken from two simulation runs 
in this paper. The effect of (3 on the flux ropes forming 
in a simple tw o-bipole configuration has previousl y been 
considered by Mackay & van Ballegooijen ( 2006a[ ). Here 
coronal flux ropes formed not only within the two bipoles 



netic measurements (Pevtsov et al.||1995 


Abramenko et 


al. 


1996 


Hagino & Sakurai 


2004J, but also by proxy 



75bservations suc h as H a image s of active regio n struc- 



ture flHal e||1927| |Balasubrama niam et al. 2004), X-ray 



sigmoids (Canfleld et al. 2007) 
measurements (Smith & Bicb< 



and in situ heliospheric 



1993). In accordance 



with this pattern we choose opposite signs of (3 in each 
hemisphere. 

To determine the effect of the (3 parameter on the mag- 



G 



Yeates et al. 



netic flux ropes forming in the global simulation, we con- 
sider in this paper a number of simulation runs, as sum- 
marized in Table [T] These include runs where the emerg- 
ing bipole twist takes the observed majority sign in each 
hemisphere (A2, A4, and A6), runs where it takes the op- 
posite sign (Am2, Am4, and Am6), and a run where the 
bipoles emerge untwisted (AO). For comparison, we also 
include a run with no emerging bipoles (AN) , a run with 
lower coronal diffusivity (D4), and a run with a lower 
outflow speed at the upper boundary (V4). 



3. FLUX ROPE FORMATION AND EVOLUTION 

In this section we outline the basic physical processes 
of magnetic flux rope formation and ejection that are 
captured in our model, illustrated with an individual flux 
rope from one of the global simulations. Overall statistics 
of flux ropes in the global simulations are presented in 
Section |5] 

Under assumptions 1 to 4 in Section \2.1\ a robust fea- 
ture of the coronal evolution is the quasi-static accumula- 
tion above photospheric polarity inversion lines (PILs) of 
axial magnetic field (i.e., the component parallel to the 
PIL). This has been demonstrated in numerical simula- 



such a field, 



tions (jya^ Ballcgooijcn fc Martens|1989||van Ballegooijen 
et al.|lM5[TLinker et al.|2003p . Field "Tines naturally con 



verge towards PILs because of surface diffusion of their 
photospheric footpoints. Although there is c ancellation 
of oppos ite polarity magnetic flux at the PIL (Martin et 
al.|1985| , any axial component of magnetic flux must re- 
main m the corona, because the photosphere p resents a 
barrier to the submergenc e of horizontal field ( van Bal- 
legooijen fc Mackay|[2007l ). The origin of the axial com- 
ponent may be simple shearing by differential rotation, 
but it may also arise from the morphology of emerging 
active regions, particularly if they contain non-zero helic- 
ity. Paper III showed that the model used in the present 
paper is able, over time, to form the hemispheric pattern 
of axial magnetic fields in filaments ( Martin et al.||1994 1 , 
precisely because it allows for these different sources of 
axial field. More generally, the formation of localized flux 
ropes in the corona is suggested to be a natural result of 
the approximate conservation of total m agnetic helicity 
durin g relaxation to a lower energy state ( Gibson & Fan 
20061 



Figure [3] illustrates the formation and ejection of a flux 
rope from one of our global simulations (run A2). Here, 
the flux rope forms above the internal PIL of a bipo- 
lar region. Panels (a) and (b) show the formation of a 
flux rope as the footpoints of sheared loops cancel at the 
PIL. The basic mech anism for formation of a twisted flux 
rope is described by Ivan Ballegooijen fc Martens] ( |1989[ |. 
Diffusivity in the corona allows two loops (such as those 
labelled X and Y in Figure |3]d) to reconnect, leading to 
the formation of longer, twisted, axial field lines such as 
that labelled Z (formed by an earlier cancellation). In 
this example, an axial magnetic component was present 
initially because the bipole was inserted with non-zero 
twist, although it has been enhanced by differential ro- 
tation over 22 days. 

In an unbounded force-free field o utside a spher e (such 
as the solar photosphere r = Rq), | Zhang et al.| (2006) 
suggest that there is an upper limit to the helicity that 
may accumulate. This arises from the virial theorem for 



r>R ( . 



2^o ' 



dV 



1 

2^o 



[Bi-Bi 



Bl)dQ, (5) 



r=R c 



which places an upper limit on the magnitude of the 
horizontal magnetic field components Bg, B^ on the sur- 
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face r = Rq, and thus on the helicity. 



( 2006 1 verify this conjecture for a particular sequence of 



axisymmctric solutions. It suggests that a flux rope can- 
not remain in force-free equilibrium once its axial field 
becomes too strong. Indeed, the sudden loss of equilib- 
rium as a flux rope gradually increases in size is a robust 
feature demonstrated both in sim ple analytical "c atas- 

and 



trophe models" (see the review by Lin et al 



in full M HD simulations (|Amari et al 
[all [20031) 



2003) 



2000 



Linker et 

Such sudden losses of equilibrium occur in 
our model when too much axial flux accumulates in flux 
ropes over the quasi-static evolution. A key aim of this 
paper is to determine what proportion of CMEs could be 
accounted for by this quasi-static model, and what pro- 
portion require some other initiation mechanism (Section 

Loss of equilibrium of a flux rope as the axial field be- 
comes too strong has previously been demonstrated in 
a two-bipole con figuration using the present magn eto- 
frictional model (Mackay & van Ballegooijen] [2006a). It 
was found that, after equilibrium is lost, the flux rope is 
driven upward by reconnection underneath, and expelled 
through the upper boundary of the computational do- 
main. Due to its simplified nature, the model does not 
follow the true dynamics upon eruption; this would re- 
quire full MHD simulations that include the inertia of the 
coronal plasma. However, we believe that the model de- 
scribes the build-up to eruption accurately. The ejection 
removes strongly-sheared field lines and magnetic helic- 
ity through the upper boundary, leaving only a small 
residual helicity at low heights. After the ejection, a 
new flux rope begins to form by flux cancellation at the 
same location, and this cycle continues. Note that the 
eruption is localized and regions away from the site of 
the eruption are unaffected and remain in equilibrium. 
Figure p3|c) shows our example flux rope just after equi- 
librium nas been lost and the axis is accelerating up- 
ward. Many of the flux rope field lines have already been 
opened (shaded in lighter grey) . Figure [3^d) shows the 
configuration after the flux rope has been ejected from 
the computational domain — note the rem aining sheared 



(but shorter) field lin es low in the corona. Mackay & van 
(j2006a| found that the frequency of lift-offs 
ile configuration depended both on the hc- 



Ball egooijen 

in their simple configuration depended 
licity of the bipolar regions and on the strength of the 
overlying field. Thus in our global simulations we expect 
the behavior of flux ropes to vary both due to the value 
of bipole twist /3, and due to the local background mag- 
netic field at each location. We now consider the global 
simulations. 

4. AUTOMATED DETECTION METHODS 

In this section, we des crib e automated techniqu es to 
detect flux ropes (Section 4.1 ) and ejections (Section 4.2 1 
in sequences of daily magnetic field snapshots from the 
global simulations. Global statistics derived with these 
techniques are then presented in Section [5] 
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Fig. 3. — Example flux rope forming in a bipolar region in simulation run A2. Top row shows the magnetic field on days (a) 191, (b) 
213, (c) 224, and (d) 226. Bottom row (panels e to h) shows vertical components of the magnetic pressure force (solid lines) and magnetic 
tension force (dashed lines) on the same days, as a function of radius above the center of the bipolar region (in arbitrary units). In panels 
(a) to (d) grey shading and thin contours show the radial magnetic field Bq t on the photosphere, and thick lines are selected coronal field 
lines (lighter grey denotes open field lines). In panels (f) and (g) the dashed vertical line indicates the height of the flux rope axis. 



4.1. Detection of Flux Ropes 

The basic idea is to search for the distinctive Lorentz 
force structure of a flux rope. This is illustrated in Fig- 
ures |3je) to (h) for the example given in Section [3] The 
Lorentz force may be decomposed as 



the magnetic pressure and tension forces are computed 
from the magnetic field B . We then require that 



Jo x B 



1 

Ho 



(B ■ V) B - V 



2/i 



(6) 



P z {r i+1 ,6 i ,(j) i )>Q, 
T z (ri-i,9i,<f>i)>0, 
T z (r i+ i,^,^)<0 



(7) 
(8) 
(9) 
(10) 



where the first term on the right-hand side is the mag- 
netic tension T and the second is the magnetic pressure 
gradient P. Figures [3je) to (h) show the vertical com- 
ponents of each of these terms T z (dashed line) and P z 
(solid line), as a function of height above the centre of 
the PIL. Notice that the two terms are approximately 
equal and opposite, consistent with a force-free equilib- 
rium with vanishing Lorentz force. The distinctive flux 
rope structure consists of a sign reversal in both P z and 
T z at the height of the flux rope axis, denoted by the ver- 
tical dashed line in Figures [3[f) and (g). In a flux rope, 
the helical structure is such that the tension force acts 
inward toward the r ope axis, but the magnetic pressure 



force acts outward (TMackay & van Ballegooijen||2006a 
see their Fig. 13) 

The automated procedure for detecting flux ropes then 
consists of two stages: 

1 . Point testing: Test individual points on the numer- 
ical grid for the characteristic flux rope structure of 
sign changes in the vertical magnetic pressure and 
tension forces. 



Clustering: Use a hierarchical clustering algorithm 
(Johnson 1967) to associate neighbouring points 
which form part of the same flux rope structure. 



for the point (r^, 9i, <pi) to be selected. In practice, we do 
not need to test all points, but rather test points only up 
to height r = 1.44i? Q , and use a coarser "testing grid" 
of (21, 93, 120) points in the (r, 9, (j>) directions. In the 
(j) direction at the equator this corresponds to one third 
of the computational grid resolution. The testing grid 
is chosen such that each point represents an equal 3D 
volume, by taking equal steps in cos 9, 4>, and r 3 . Using 
a coarser grid reduces computational effort and does not 
affect the results since we are only interested in well- 
resolved flux rope structures. In order to concentrate 
the sample on twisted flux ropes, we have implemented 
the fifth condition of a sufficiently strong parallel current 
at each point tested, requiring that 



|jo-Bo| 

B 2 

D 



> a* 



(11) 



First consider stage 1. At each point (r^, 9t 1 tfii) in the 
computational grid the vertical components P z and T z of 



with the threshold a* = 0.7 x 10 _8 m _1 determined by 
experiment. By way of example, Figure |4^a) shows the 
flux rope points identified by this first stage on day 202 of 
run A2, mid- way through the simulation. The points are 
projected on a plot of B$ r on the photosphere, showing 
the PIL dividing regions of positive (white) and negative 
(grey) magnetic polarity. The identified points tend to 
lie above PILs, as expected from our basic theory of flux 
rope formation (Section [3]) . Furthermore, the points are 
mostly grouped into larger structures. Automated detec- 
tion of these groupings forms stage 2 of the procedure. 
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structure. 



4.2. Detection of Ejections 
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Fig. 4. — Automated detection of magnetic flux ropes in the 
simulated magnetic field, for day 202 of run A2. The first two 
panels show (a) selected points and (b) clusters, projected on the 
normal magnetic field Bq t in the photosphere (white shading for 
positive, grey for negative). Panel (c) shows coronal magnetic field 
lines traced from the selected flux rope points (with lighter grey 
denoting open field lines), viewed from 180° longitude and 20° 
latitude. 



The basic clustering idea is simple. Starting with each 
point as an individual entity, the two closest points are 
grouped together, and the process repeated until the 
shortest inter-group distance (in 3D space) is above some 
threshold. For this threshold we choose 5Ai? , where A 
is the heliographic angle in radians of a computational 
grid cell at the equator. After running this clustering al- 
gorithm, any groups with fewer than 8 points on the test- 
ing grid are removed. Again this value was determined by 
visual inspection of the selected structures. The results 
after clustering for day 202 of run A2 are shown in Figure 
EFb). Each individual group of points, or "flux rope", is 
identified by color and a number. The actual magnetic 
field structures corresponding to these flux ropes are il- 
lustrated in Figure EFc) , which shows 3D magnetic field 
lines traced from trie selected flux rope points in each 



An immediate problem is how to define a flux rope 
ejection within our simulations. We adopt the practi- 
cal definition of a large enough radial velocity in the 
magneto-frictional code. Thus we include both losses 
of equilibrium following gradual build-up of axial field, 
and sudden rises caused by nearby bipole emergence. We 
also include partial lift-offs where only one end of a flux 
rope opens up. Typically the other end is held down 
by overlying field from nearby regions. The majority of 
ejections remove the main, twisted, part of the flux rope 
through the top boundary of the computational box, usu- 
ally within several days of the onset of rapid acceleration. 

Our automated procedure to find ejections is straight- 
forward: 

1. Record the radial velocity u ur in the magneto- 
frictional code at each of the selected flux rope 
points over the simulation, and select those points 
where v 0r > 0.5 km s -1 . 

2. Cluster these points into separate ejection events. 

Although a similar algorithm is used, this clustering pro- 
cedure differs from that used to define flux ropes on any 
particular day. Firstly, the clustering is carried out both 
in space and time, as ejections can take up to a few days 
to leave the computational domain. Secondly, for the 
distance measure we consider only the heliocentric angle 
between points, neglecting any radial distance. This al- 
lows for radial movement between daily snapshots. We 
require a minimum inter-group angle of 4A between clus- 
ters, and a minimum separation in time of 4 days. After 
clustering, all groups of points with fewer than 4 points 
are discounted. These optimal parameter values have 
been determined by testing the detection procedure over 
a 13-day test period from day 202 to day 214 in run A2. 
Locations and days of detected ejections in this test pe- 
riod were compared to those found in a careful manual 
study of the magnetic field structure. The performance 
of the automated detection algorithm in this period (with 
the optimal choice of parameters found) is demonstrated 
by Table [2] and also illustrated in Figure [5 a), which 
shows all nux rope points during the period superim- 
posed on the normal photospheric magnetic field for day 
208. Panels (b), (c), and (d) correspond to oth er simula- 
tion runs and will be discussed in Section 15.21 The first 
column of Table |2]lists the days of significant vertical mo- 
tion of the flux rope in each actual ejection, as detected 
manually. The approximate locations in latitude and 
longitude arc listed in the second column, and plotted 
as black triangles in Figure [5j a). The third and fourth 
columns of Table [2] give the results of the automated de- 
tection algorithm (with the optimal parameters). The 
clusters of points selected in this procedure are colored 
red in Figure |5|a) . Note that only points falling within 
the 13-day period are shown in the figure, even if the 
cluster overlaps with the period. In Table [2j the num- 
ber of clustered points refers to the total number in that 
cluster, whether in the 13-day period or not. The code 
detects 16 of the 19 ejections (84%), and the number of 
clustered points (fourth column) gives an estimate of the 
size of each event. Of the three missed ejections, two are 
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TABLE 2 

Flux rope ejections in run A2 between day 202 and day 214. 



Actual Days LjcitituQG o 


T '+ A [A \a 

z Longitude ^dog.J 


Detected Days 


Clustered Joints 


198—203 


45, 185 


195—202 (197) 


60 


202 


15, 245 


195—202 (197) 


135 


202—203 


10, 315 


201—203 (202) 


16 


203—204 


5, 230 


200—205 (205) 


28 


204—206 


65, 200 


202—204 (204) 


46 


205 


-40, 260 


not detected 




205—206 


-20, 195 


205 


19 


206 


60, 175 


206—207 (206) 


16 


206—207 


-35, 270 


206—207 (206) 


10 


207—208 


-20, 175 


207—209 (208) 


10 


208 


25, 55 


208 


7 


208—212 


-60, 5 


not detected 




209 


35, 105 


207—209 (209) 


28 


209—211 


40, 340 


208—212 (211) 


204 


211—212 


45, 170 


not detected 




212—213 


-30, 125 


212 


24 


214 


-20, 285 


215—216 (215) 


99 


214+ 


-20, 310 


213—217 (217) 


63 


214+ 


55, 15 


210—215 (214) 


32 


False detections: 




201—204 (201) 


17 






205—206 (205) 


10 






207—208 (207) 


10 






212 


■1 



a Approximate location of manually-detected ejection. 15 Day with most selected points is 
given in parentheses. c Total number of selected points in detected ejection cluster over all 
days. 
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Fig. 5. — Locations of flux rope ejections during the 13-day period from day 202 to day 214, used to test the automated detection 
algorithm. All flux rope points during the period are shown for runs (a) A2, (b) Am2, (c) A4, and (d) D4. Points are colored red if 
any point at that (8, <j>) location was selected as part of an ejection by the algorithm, blue otherwise. The points are projected on to 
the normal photosphcric magnetic field Bo r from day 208 (white for positive and grey for negative). In panel (a), black triangles show 
manually-determined locations of ejections for run A2. 
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small, and the third (located at -60, 5) is quite large but 
rises slowly and gradually, with no sudden loss of equilib- 
rium. We note that there are also four "false positives," 
where spurious ejections are detected. However, two of 
these are small "flux rope" structures associated with 
temporary coronal reconfiguration following new bipole 
emergence, and the other two were actually considered 
on manual inspection to be part of other ejections. We 
conclude that this procedure performs well and the num- 
ber of ejections detected is accurate to within ±15%. 

5. RESULTS 

We have applied the automated detection algorithms 
described in Section [4] to each of the 177-day simulation 
runs in Table [TJ In this section we discuss the results, 
concentrating first on over all properties of the flux ropes 
in each run (Section 5.1|, a nd th en on the number of 
ejections in each run (Section 5.2 1. 

To give an overall feeling for the evolution, Figure [6] 
shows a sequence of snapshots of the flux ropes in sim- 
ulation run A2 at 27-day inter vals , as detected by the 
automated procedure in Section 4.1 Note that the num- 
ber and size of flux ropes increase early in the simulation 
as the magnetic field diverges from the initial potential 
field, before the numbers begin to saturate. Their loca- 
tions evolve as new bipoles emerge and the coronal field 
is reconfigured by surface motions. 

5.1. Dependence of Flux Rope Properties on 
Simulation Parameters 

Figure [7] shows how the properties of flux ropes evolve 
over time in each simulation run. These properties in- 
clude the total number of flux rope points (Figure [7b,) , 
the mean latitude of flux rope points (Figure 17b) , the 
number of flux ropes found by the clustering stage (Fig- 
ure 7;), and the mean size of each of these clusters (Fig- 
ure 71). 

A general feature in all of these properties is an ini- 
tial increase, followed in most cases by a saturation after 
50 or more days. This illustrates the gradual build-up 
of hclicity in the corona as the simulation evolves away 
from the initial potential field. Because run AN (with no 
bipole emergence) also shows an initial increase in the 
number of flux rope points (Figure , we see that sur- 
face shearing alone creates flux ropes over this timescale. 
However, the faster increase for the other runs shows that 
emerging helicity speeds up the formation of flux ropes. 
Indeed, Figure mb) shows that more flux ropes form at 
lower latitudes m the runs with bipole emergence, sup- 
porting this point. The emergence of new flux also acts 
to sustain the coronal magnetic field against the diffu- 
sive decay which would otherwise occur over hundreds 
of days. This decay is evident in run AN as reformation 
of ejected flux ropes and production of new flux ropes 
declines after day 210. 

Notice that all runs show a temporary decrease in the 
number of flux rope points and in their mean latitude 
around day 210. This is caused by the ejection of a pair 
of high- latitude flux ropes, illustrated in a time sequence 
in Figure [8] (for run A2) . The two ejections are labelled 
El and E2, and take place from days 199 and 204 re- 
spectively. Ejection E2 is also seen in Figure pna), but 
El occurred before the period shown there. A similar 



high-latitude ejection causes a dip in the number of flux 
rope points around day 260 for run AN. 

5.1.1. Effect of Emerging Bipole Twist 

Consider the effect of the emerging bipole twist (3, as 
shown by colors in Figure [7] Firstly, the number of flux 
ropes present at any time (Figure [7b) saturates at around 
38 in all runs, although each shows considerable fluctu- 
ation between about 28 and 48. There is some tendency 
toward more flux ropes in runs A6 and Am6, but this is 
not clear at all times. Run AO, where bipoles emerge un- 
twisted, still has a comparable number of flux ropes, thus 
suggesting that the maintenance of photospheric flux is 
sufficient to produce flux ropes, without requiring the 
emergence of already-twisted fields. 

There are, however, differences in the mean latitude 
and mean size of flux ropes in runs with different (3. Over 
the full simulation, the mean latitude is lower for runs 
with more emerging bipole twist (28.8° for run A6 and 
27.3° for run Am6 versus 31.1° for run AO). This sug- 
gests that greater emerging helicity accelerates flux rope 
formation at lower latitudes, although it doesn't signifi- 
cantly affect the number of locations at which flux ropes 
form. 

Figure [7 d) shows that, later in the simulation, the flux 
ropes arelarger in the runs with majority hemispheric 
sign of bipole twist (A2, A4, and A6) as compared to the 
runs with opposite sign of twist (Am2, Am4, Am6). A 
more detailed analysis shows that the flux ropes that are 
larger in the first three runs are located mainly at mid- 
latitudes (40° to 60°), rather than at high latitudes or 
among newly-emerged active regions. The explanation 
is suggested by our analysis of sheared magnetic fields 
in filaments (Paper III). In particular, consider the PIL 
internal to a single bipole. Initially, the same amount of 
shear will emerge whether (3 — —0.6 or j3 — 0.6, although 
with opposite sign. Hence there is no difference in flux 
rope sizes between runs Am6 and A6 at active latitudes. 
However, differential rotation always acts to produce the 
observed majority chirality (direction of shear) over such 
north-south oriented PILs. In run A6, the emerged shear 
already has the majority direction, but in run Am6 the 
emerged shear has the opposite direction, so the shear 
will first be reduced before the direction is reversed. Thus 
at mid-latitudes, i.e., among older regions, the amount 
of shear and hence the size of the flux ropes will be lower 
in run Am6 than in run A6. There is no difference at 
high latitudes because, in the 177-day simulations, these 
are not influenced by the emerging regions. 

In summary, the main effects of emerging bipole twist 
are to speed production of flux ropes at active latitudes 
for greater magnitude of (3, and to produce larger flux 
ropes at mid-latitudes for the correct sign of j3. The 
number of flux ropes present at any time is not affected 
in a systematic way. 

5.1.2. Effect of Radial Outflow 

In run V4, shown by the thin solid black line in Figure 
[7J we halve the magnitude of the radial outflow speed to 
50 km s -1 . We find no significant effect on the overall flux 
rope properties. This is because the outflow speed has 
no direct influence on the formation of flux ropes, or on 
the initial stages of the ejection. It acts only at the top 
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Fig. 6. — Evolution of flux ropes during run A2, showing results of the automated procedure on six representative days. As in Figure |4[b), 
detected ropes are shown by colored points, projected on the normal magnetic field Bq t in the photosphere (white shading for positive, 
grey for negative). 



of the computational domain, so influences only the (un- 
physical) evolution after equilibrium has been lost. To 
minimize dependence on this unphysical evolution after 
loss of equilibrium, our automated detection algorithm 
considers only flux ropes below r = 1.44i? Q . 

5.1.3. Effect of Coronal Diffusion 

The dashed black line in Figure [7] shows run D4, where 
the coronal diffusivity 770 nas only half of its usual value. 
This results in 20% more flux rope points being found 
(Figure [7k) , and in 20% larger flux ropes (|7ji) . This is 
because the coronal diffusion acts to dissipate the con- 
centrated currents in twisted flux rope structures. Lower 
diffusion allows more highly twisted structures to form. 
However, the coronal diffusion is not the primary source 
of axial magnetic flux above PILs, thus does not deter- 
mine the locations at which flux ropes form. Figure [7jc) 
shows that the mean latitude of flux ropes is not affected 
by the coronal diffusivity. A natural result of the larger, 
more twisted flux ropes formed in run D4 is to produce 
more ejections. We consider the number of ejections pro- 



duced in different runs next. 

5.2. Rate of Flux Rope Ejections 

Using the automated procedure in Section [4. 2[ we find 
the number of ejections per day to have a ramp-up phase 
of about 40 days, before fluctuating considerably over 
the rest of the simulation, as well as varying between 
runs. Figure [9] shows the total number of ejections in 
each run, with error bars showing the estimat ed e rrors 
in the automated detection procedure (Section |4.2| . We 
immediately see the importance of emerging bipolcs for 
flux rope ejections. In run AN there are only 23 ±4 ejec- 
tions, as compared to 108 ± 16 in run AO, and more for 
the runs with emerging bipole twist — up to 202 ± 30 for 
run A6. Partly, this reflects the increased number of flux 
ropes in runs with emerging bipoles, and partly it re- 
flects the influence of strong emerging magnetic fields in 
reconfiguring the coronal field. The increasing number of 
ejections with greater emerging bipole twist reflects the 
quicker formation of flux ropes with strong axial fields, 
which then lose equilibrium. For example, compare pan- 
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TABLE 3 
Rates of simulated ejections and 
observed cmes. 



Fig. 7. — Dependence of flux rope statistics on simulation pa- 
rameters: (a) total number of points selected (on testing grid) after 
clustering, (b) mean latitude of selected points, (c) number of flux 
ropes present, and (d) mean number of points (on testing grid) per 
flux rope. Colored lines show runs with different emerging bipole 
twist /3 (see legend in panel a). Thick black line shows run AN, thin 
solid black line shows run V4, and thin dashed black line shows run 
D4. In panel (b) vertical bars show la for run A4. 



els (a) and (c) in Figure [5] which show flux rope points 
and radial velocities for runs A2 and A4. The flux rope 
labelled M in Figure p[c), internal to a newly-emerged 
bipole, both forms andis ejected in run A4. However, in 
run A2 there is not sufficient axial magnetic field at this 
location for a flux rope to be detected at this time. 

There is also a tendency for the runs with the majority 
hemispheric sign of bipole twist (A2, A4, and A6) to have 
more ejections than the runs with the same magnitude 
but opposite sign of twist. This may be related to the 
larger size of mid-latitude flux ropes in the first set of 



Simulation Run 


Ejections per Day a 


AN 


0.09 ± 0.01 


Am6 


1.05 ± 0.16 


Am4 


0.88 ±0.13 


Am2 


0.61 ± 0.09 


AO 


0.67 ± 0.10 


A2 


1.01 ± 0.15 


A4 


1.20 ±0.18 


A6 


1.28 ±0.19 


Dl 


1.60 ±0.24 


V4 


1.16 ±0.17 


LASCO events: 




All 


3.59 


"Poor events" removed 


2.25 



a Between day 183 (1999 July 2) and day 283 
(1999 October 10). 



runs, as described in Section 5.1 This is illustrated by 
comparing run Am2 with run A2 in Figure [5] (panels 
b and a respectively) . The black boxes in Figure [5jb) 
show decaying bipolar regions where no flux ropes are 
detected in run Am2. In run A2, both of these locations 
have already formed strong flux ropes which are ejected 
during this particular correlation period. 

From run V4 — shown by a square in Figure [9] — we see 
that halving the outflow speed has no significant influ- 
ence on the number of ejections (giving 178±27). So, just 
as the outflow has no major influence on the formation 
of flux ropes, it has no major influence on their loss of 
equilibrium. As before, this is because it acts only near 
the upper boundary, which a flux rope reaches only after 
equilibrium is lost. In contrast, halving the coronal dif- 
fusivity 770 increases the number of ejections to 233 ± 35 
in run D4 (triangle in Figure |9| . An example of an ad- 
ditional ejection not occurring m run A4 is labelled N in 
Figure p[d). The higher number of ejections is explained 
by the larger, more twisted flux ropes which are able to 
form in this case. 

Finally, how does the number of ejections produced 
in our simulations compare with the number of CMEs 
observed on the real Sun duri ng this period? The CD AW 
catalog ( |Yashiro et al. 20041 is the standard manually- 
compiled list of CMEs observed by the SOHO/LASCO 
coronagraphs (Brueckner et al. 1995 1. To avoid the initial 



ramp-up phase in the simulation, we compare the rates 
of ejections between day 183 (1999 July 2) and day 283 
(1999 October 10). Table [3] shows the number of flux 
rope ejections per day for each simulation run over this 
period, in addition to the observed rate from the CD AW 
catalog. Two observed rates are given: the first includes 
all events in the catalog, and the second omits events 
labelled "poor" by the LASCO operator. We believe the 
second rate to be more appropriate for this comparison, 
as our global simulations relate to large-scale flux rope 
events, which are unlikely to be labelled "poor". This 
rate will be a lower estimate for two reasons: firstly not 
all far-side events are observed, and, secondly, there are 
several gaps in the instrument coverage over the period. 
Using this rough estimate, our simulations produce flux 
rope ejections at about 50% of the observed CME rate 
(after the initial ramp-up phase). 
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Fig. 8.— Sequence of two high-latitude ejections (El and E2) in run A2, seen on days (a) 197, (b) 200, (c) 202, (d) 204, (e) 205, and 
(f) 206. Grey shading and thin contours show radial magnetic field Bq t on photosphere (white/solid contours for positive, grey/dashed 
contours for negative). Thick lines show selected coronal magnetic field lines traced from flux rope points, with lighter grey denoting open 
field lines. 




-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 
Emerging Bipolc Twist (S. Hemisphere) 

Fig. 9. — Total number of ejections in different simulation runs, 
as detected by the automated procedure. Asterisks joined by solid 
line show runs Am6, Am4, Am2, AO, A2, A4, and A6. Dashed line 
shows run AN (where emerging bipolc twist is irrelevant), triangle 
run D4, and square run V4. Error bars show estimated error in 
automated detection algorithm. 

6. CONCLUSION 

We have studied the formation and ejection of mag- 
netic flux ropes in a simplified model of the coronal mag- 



netic field evolution, to begin to address the question of 
where and when CMEs occur in the global context. Loss 
of equilibrium of magnetic flux ropes in the low corona 
is a leading model for the initiation of CMEs, and the 
model described in this paper is, in essence, an extension 
of previous 2D catastrophe models of a single eruption to 
the 3D global corona. Ultimately, we aim to determine 
what proportion of observed CMEs may be explained by 
the loss of equilibrium mechanism. Using automated de- 
tection algorithms, we have tracked the formation, loss of 
equilibrium, and ejection of flux rope structures forming 
at many locations in the simulated corona, which evolves 
continually in response to the emergence of twisted bipo- 
lar active regions, and to large-scale motions on the pho- 
tospheric boundary. 

In this paper we consider the effect of key simulation 
parameters on the flux ropes formed and on the rate of 
flux rope ejections. We draw the following main conclu- 
sions: 

1. The number of flux ropes present at any one time 
fluctuates between about 28 and 48, with no sys- 
tematic dependence on the hclicity of emerging 
bipolcs. If no new bipoles are emerged, the surface 
flux decays and the number of flux ropes decreases 
due to ejections. 

2. The magnitude of emerging bipolc hclicity has no 
major effect on the number of flux ropes present at 
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any one time, but greater emerging helicity leads 
to more flux rope ejections. 

3. The sign of emerging helicity also has an effect. If 
active regions emerge with the (observed) minority 
sign of twist in each hemisphere, then smaller flux 
ropes are created at mid-latitudes, and there are 
fewer ejections. This is because the surface shear- 
ing has first to reverse the direction of the sheared 
field that emerged with the minority sign of helic- 
ity, before forming a flux rope with the majority 
helicity. 

4. The results are not sensitive to the upper bound- 
ary condition (radial outflow), but do depend on 
the turbulent diffusivity in the corona. A lower 
diffusivity leads to larger flux ropes and more ejec- 
tions. 

5. The rate of flux rope ejections in the model is 
roughly 50% of the observed LASCO CME rate, de- 
pending on the choice of parameters in the model. 

We have shown that the rate of flux rope ejections 
varies from 0.67 ± 0.10 per day in run AO (when bipolcs 
emerge untwisted) to 1.28 ± 0.19 in run A6 (the great- 
est amount of emerging helicity considered). Since we 
do not, at present, have reliable measurements of the 
magnetic helicity in all 119 active regions that emerged 
during the simulated period, we cannot predict this ejec- 
tion rate precisely. However, constraints from filament 
chirality observations (Paper II) suggest that simulation 
runs A4 and A6 are likely to be most realistic. Hence 
our conclusion that the quasi-static model can produce 
about 50% of observed CMEs. Since the sign and magni- 
tude of helicity emerging in active regions play an impor- 
tant role, more realistic future models must ultimately 
incorporate observations of helicity in individual emerg- 
ing regions, which may come from the forthcoming Solar 
Dynamics Observatory (SDO) mission. Note that the 
rate of ejections in the model may be increased some- 
what by reducing the turbulent diffusivity in the corona, 
although too small a value (as in run D4) produces flux 
ropes which are more highly twisted than typically ob- 
served. Future simulations will consider the effect of a 
higher-order form of "hyperdiffusion," which has been 
suggested to be more appropriate for the solar corona 
because it conser ves magnetic helicity (van Ballegooijen 
fc Cranmer|[2008 l 



the amount of axial magnetic flux removed in ejections, 
and that in reality more axial flux is left behind, leading 
to another ejection from the same location shortly after. 
However, the axial flux rebuilds over a period of tens 
of days following an ejection, so this quasi-static model 
may never be able to produce a rapid succession of flares 
and eruptions from the same active region over minutes 
to hours, as is sometimed observed (Gopalswamy et al. 



2006 1. It seems likely, therefore, that up to 507o of CMEs 
require either more complex magnetic field structures 
than allowed for in our simple model, or some physical 
mechanism other than loss of equilibrium of flux ropes 
formed quasi-statically by flux cancellation. One possi- 
bility is the dynamic emergence of a lready twiste d flux 
ropes, and subsequent activity (e.g., Tanaka 1991|. Al- 



ternatively, in full MHD , eruptions may be triggered by 
dTokman & Bellan||2002| f 



& Kliem|2003 


Kusano et al.|2004||Kliem & T6r6k|2006p 


Finally, in tl 


lis paper we have concentrated on the ejec 



Why does the model produce only 50% of the observed 
CME rate? It is possible that the model over-estimates 



tion of flux ropes within the model, making only a rough 
comparison with the rate of observed CMEs. A future 
paper will describe a more detailed comparison with the 
source regions of observed CMEs, in order to determine 
whether the simplified global model can help to explain 
or predict the initiation locations of these events. In 
addition, it will discuss how meaningful present compar- 
isons to observations are, and what extra input is re- 
quired from observations so that more meaningful com- 
parisons may be made. Only through such studies will 
we gain the insight of what is needed to predict CMEs 
in future. 
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